!*******************************************************************************************************************
! This file defines fit function objects. It will fit rational function
! approximations and chebyshev polynomial approximations.
!
! The idea is to declare a fit type:
!
! type(tRFit) :: fit --- for rational function
! type(tCFit) :: fit --- for chebyshev polynomial
!
! then use new to do the fit:
!
! call new(fit, ll, ul, [inverse], [xLog], [yLog], [param], [universe], [lps])
!
! with user supplied lower and upper limits for the fit. There are other parameters you can change through the new
! method, like orders of the polynomials and such. You can find the functions below (new_rFit and new_cFit).
!
! [inverse] --- set to .true. to also fit the inverse of the function.
! [xLog], [yLog] --- .true. changes the variable to log units. This happens internally, so you don't have to add log or exp
!                     to your code.
! [param] (real(D)), [universe], [lps] --- pass these types into your fit function.
!
!*** ---> evaluate(fit, x)
!
! A good example is everything in base/nu.F90.
!
!******************************************************************************************************************
module mFit

        use mPrecision
        use mObject
        use mUniverse
        use mLPS

        implicit none

        ! Defaults for rational function fit
        integer, parameter :: DEFAULT_NUM_ORDER = 5, &
                              DEFAULT_DEN_ORDER = 5
        ! # of Chebyshev coeffs for fit
        integer, parameter :: NPFAC = 8 ! numerical recipes suggests 8

        interface new
                module procedure new_rFit
                module procedure new_cFit
        end interface

        interface param
                module procedure param_fitParams
        end interface param

        interface getLPS
                module procedure getLPS_fitParams
        end interface getLPS

        interface getUniverse
                module procedure getUniverse_fitParams
        end interface getUniverse

        interface fit
                module procedure fitRationalFunction
                module procedure fitChebyshev
        end interface fit

        interface writeCoeffs
                module procedure writeRFitCoeffs
                module procedure writeCFitCoeffs
        end interface writeCoeffs

        interface switchPointers
                module procedure switch1DPointers
                module procedure switch2DPointers
        end interface switchPointers

        interface evaluate
                module procedure evalRationalFunction
                module procedure evalChebyshevFunction
        end interface evaluate

        interface changeX
                module procedure rFitChangeX
                module procedure cFitChangeX
        end interface changeX

        interface changeY
                module procedure rFitChangeY
        end interface changeY

        interface fitLL
                module procedure cFitLL
                module procedure rFitLL
        end interface fitLL

        interface fitUL
                module procedure cFitUL
                module procedure rFitUL
        end interface fitUL

        interface evaluateInverse
                module procedure fitEvaluateInverse
        end interface evaluateInverse

        interface fitDone
                module procedure rFitDone
                module procedure cFitDone
        end interface fitDone

        private
        public :: new
        public :: param, getLPS, getUniverse, getLPSUniverse
        public :: writeOrders
        public :: writeRFitCoeffs,writeCFitCoeffs
        public :: fit,evalRationalFunction,fitEvaluateInverse,evaluateInverse
        public :: fitLL,fitUL
        public :: cFitLL,rFitLL,cFitUL,rFitUL
        public :: rFitDone,cFitDone
        public :: chebyshevDerivative,evalChebyshevFunction,evaluate

contains

        pure subroutine new_rFit(rFit, ll, ul, numOrder, denOrder, &
                                 inverse, xLog, yLog, param, universe, lps)
                type(tRFit), intent(out) :: rFit
                real(D), intent(in) :: ll, ul
                integer, intent(in), optional :: numOrder, denOrder
                logical, intent(in), optional :: inverse, xLog, yLog
                real(D), intent(in), optional :: param
                type(tUniverse), intent(in), optional :: universe
                type(tLPS), intent(in), optional :: lps

                rFit%ll = ll
                rFit%ul = ul

                if( present(numOrder) ) then
                        rFit%numOrder = numOrder
                else
                        rFit%numOrder = DEFAULT_NUM_ORDER
                end if

                if( present(denOrder) ) then
                        rFit%denOrder = denOrder
                else
                        rFit%denOrder = DEFAULT_DEN_ORDER
                end if

                if( present(inverse) ) then
                        rFit%inverse = inverse
                else
                        rFit%inverse = .false.
                end if

                if( present(xLog) ) then
                        rFit%xLog = xLog
                else
                        rFit%xLog = .false.
                end if

                if( present(yLog) ) then
                        rFit%yLog = yLog
                else
                        rFit%yLog = .false.
                end if

                rFit%fitDone = .false.

                if( present(param) ) rFit%params%param = param
                if( present(universe) ) rFit%params%universe = universe
                if( present(lps) ) rFit%params%lps = lps

        end subroutine new_rFit

        real(D) pure function param_fitParams(fitParams)
                type(tFitParams), intent(in) :: fitParams
                param_fitParams = fitParams%param
        end function param_fitParams

        type(tLPS) pure function getLPS_fitParams(fitParams)
                type(tFitParams), intent(in) :: fitParams
                getLPS_fitParams = fitParams%lps
        end function getLPS_fitParams

        type(tUniverse) pure function getUniverse_fitParams(fitParams)
                type(tFitParams), intent(in) :: fitParams
                getUniverse_fitParams = fitParams%universe
        end function getUniverse_fitParams

        type(tUniverse) pure function getLPSUniverse(fitParams)
                type(tFitParams), intent(in) :: fitParams
                getLPSUniverse = getUniverse( getLPS(fitParams) )
        end function getLPSUniverse

#ifdef JUNK
        real(D) pure function param_RFit(rFit)
                type(tRFit), intent(in) :: rFit
                param_RFit = rFit%params%param
        end function param_RFit

        type(tLPS) pure function getLPS_rFit(rFit)
                type(tRFit), intent(in) :: rFit
                getLPS_rFit = rFit%params%lps
        end function getLPS_rFit
#endif JUNK

        subroutine writeRFitCoeffs(rFit)
                type(tRFit), intent(in) :: rFit
                IF(rFit%fitDone) then
                        write(0, *) rFit%coeffs
                        write(0, *)
                        IF(rFit%inverse) write(0, *) rFit%invCoeffs
                else
                        write(0,*) 'fit not done.'
                end if
        end subroutine writeRFitCoeffs

        subroutine writeOrders(rFit)
                type(tRFit), intent(in) :: rFit
                write(6, *) 'numerator order:', rFit%numOrder
                write(6, *) 'denominator order:', rFit%denOrder
        end subroutine writeOrders

        real(D) pure function rFitLL(rFit)
                type(tRFit), intent(in) :: rFit
                rFitLL = rFit%ll
        end function rFitLL

        real(D) pure function rFitUL(rFit)
                type(tRFit), intent(in) :: rFit
                rFitUL = rFit%ul
        end function rFitUL

        real(D) pure function cFitLL(cFit)
                type(tCFit), intent(in) :: cFit
                cFitLL = cFit%ll
        end function cFitLL

        real(D) pure function cFitUL(cFit)
                type(tCFit), intent(in) :: cFit
                cFitUL = cFit%ul
        end function cFitUL

        pure subroutine new_cFit(cFit, ll, ul, numCoeffs, xLog, yLog, &
                                 param, lps)

                type(tCFit), intent(out) :: cFit
                real(D), intent(in) :: ll,ul
                integer, intent(in), optional :: numCoeffs
                logical, intent(in), optional :: xLog, yLog
                real(D), intent(in), optional :: param
                type(tLPS), intent(in), optional :: lps

                cFit%ll = ll
                cFit%ul = ul
                if( present(numCoeffs) ) then
                        cFit%numCoeffs = numCoeffs
                else
                        cFit%numCoeffs = MAX_GAUSS_COEFFS
                end if
                if( present(xLog) ) then
                        cFit%xLog = xLog
                else
                        cFit%xLog = .false.
                end if
                if( present(yLog) ) then
                        cFit%yLog = yLog
                else
                        cFit%yLog = .false.
                end if
                cFit%fitDone = .false.

                if( present(param) ) cFit%params%param = param
                if( present(lps) ) cFit%params%lps = lps

        end subroutine new_cFit

#ifdef JUNK
        real(D) pure function param_cFit(fit)
                type(tCFit), intent(in) :: fit
                param_cFit = fit%params%param
        end function param_cFit

        type(tLPS) pure function getLPS_cFit(cFit)
                type(tRFit), intent(in) :: rFit
                getLPS_rFit = rFit%params%lps
        end function getLPS_rFit
#endif JUNK

        logical pure function rFitDone(rFit)
                type(tRFit), intent(in) :: rFit
                rFitDone = rFit%fitDone
        end function rFitDone

        logical pure function cFitDone(cFit)
                type(tCFit), intent(in) :: cFit
                cFitDone = cFit%fitDone
        end function cFitDone

        subroutine writeCFitCoeffs(cFit)
                type(tCFit), intent(in) :: cFit
                if(cFit%fitDone) then
                        write(0, *) cFit%coeffs
                        write(0, *)
                else
                        write(0, *) 'fit not done.'
                end if
        end subroutine writeCFitCoeffs

        !***********************************************************************
        ! Given a function func, lower and upper limits of the interval [ll,ul],
        ! this routine calculates the coefficients c_k such that
        ! func(x) ~ [Sum c_k T_(k-1)(y)]-c_1/2, where y is related to x by a
        ! change of variables. This routine should be used with a large value
        ! of n (size of the coeffs array), such as 30 or 50. The array of
        ! coeffs(:) can be truncated to a smaller value m such that c_(m+1) and
        ! subsequent elements are negligible.
        !***********************************************************************
        subroutine fitChebyshev(cFit, func)

                use mConstants

                type(tCFit), intent(inout) :: cFit
                interface
                        function func(params, x)
                                use mPrecision
                                use mObject
                                real(D) :: func
                                type(tFitParams) :: params
                                real(D), intent(in) :: x
                        end function func
                end interface
                integer :: info, k, j
                real(D) :: ll, ul, bma, bpa, floatNumCoeffs, y, xVal, fac, &
                           f(cFit%numCoeffs), sum

                ll = changeX(cFit, cFit%ll)
                ul = changeX(cFit, cFit%ul)

                bma = 0.5_D * (ul - ll)
                bpa = 0.5_D * (ul + ll)

                ! The function is evaluated at the n points required by (5.8.7)
                ! in Numerical Recipes in Fortran:
                floatNumCoeffs = real(cFit%numCoeffs, D)
                eval_func: do k = 1, cFit%numCoeffs
                        y = cos( PI * (k - 0.5_D) / floatNumCoeffs )
                        xVal = y * bma + bpa
                        if(cFit%xLog) then
                                xVal = exp(xVal)
                        end if
                        f(k) = func(cFit%params, xVal)
                        !print *,'f(',k,'; x = ',xVal,') = ',f(k)
                        !if (k==27) then
                          !print *,'xxx'
                        !endif

                        if(cFit%yLog) then
                                f(k) = log(f(k))
                        end if
                end do eval_func
                fac = 2.0_D / floatNumCoeffs
                do j = 1, cFit%numCoeffs
                        sum = 0.0_D
                        do k = 1, cFit%numCoeffs
                                sum = sum + f(k) &
                                        * cos( (PI * real( (j-1), D )) &
                                        *( (k - 0.5_D) / cFit%numCoeffs ) )
                        end do
                        cFit%coeffs(j) = fac * sum
                end do
                cFit%fitDone = .true.

        end subroutine fitChebyshev
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

        real(D) function evalChebyshevFunction(fit,x)
                type(tCFit), intent(in) :: fit
                real(D), intent(in) :: x
                integer :: j
                real(D) :: ll, ul, changedX, dHere, dd, sv, y, y2

                sanity: if(.not. fit%fitDone) then
                        write(0, *) 'ERROR: Chebyshev function not yet fit!'
                        stop
                end if sanity
                if(x < fit%ll .or. x > fit%ul) then
                        write(0, *) 'WARNING: x not in range for Chebyshev'
                        write(0, '(a, e12.5, a, e12.5, a, e12.5)') &
                                    '         Function:', fit%ll, ' <=', x, &
                                                   ' <=', fit%ul
                end if

                dHere = 0.0_D
                dd = 0.0_D

                ll = changeX(fit, fit%ll)
                ul = changeX(fit, fit%ul)
                changedX = changeX(fit, x)

                ! Change of variables:
                y = (2.0_D * changedX - ll - ul) / (ul - ll)
                y2 = 2.0_D * y
                do j = fit%numCoeffs, 2, -1 ! Clenshaw's recurrence
                        sv = dHere
                        dHere = y2 * dHere - dd + fit%coeffs(j)
                        dd = sv
                end DO
                ! The last step is different:
                evalChebyshevFunction = y * dHere - dd + 0.5_D * fit%coeffs(1)
                !print *,'Eval(',x,') = ',evalChebyshevFunction
                !print *,y,dHere,dd,fit%coeffs(1)
                !print *,'exp = ',exp(evalChebyshevFunction)
                if(fit%yLog) evalChebyshevFunction = exp(evalChebyshevFunction)
                
        end function evalChebyshevFunction
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

        type(tCFit) function chebyshevDerivative(fit)
                type(tCFit), intent(in) :: fit
                integer :: info
                integer :: j
                real(D) :: ll, ul, con

                sanity: if(.not. fit%fitDone) then
                        write(0, *) 'ERROR: Tried to differentiate an undone'
                        write(0, *) '       Chebyshev fit.'
                        stop
                end if sanity

                call new(chebyshevDerivative, fit%ll, fit%ul, fit%numCoeffs, &
                         fit%xLog, fit%yLog)

                chebyshevDerivative%coeffs(chebyshevDerivative%numCoeffs) = &
                        0.0_D
                chebyshevDerivative%coeffs(chebyshevDerivative%numCoeffs - 1) &
                        = 2.0_D * (chebyshevDerivative%numCoeffs - 1) &
                                * fit%coeffs(chebyshevDerivative%numCoeffs)
                if(chebyshevDerivative%numCoeffs > 3) then
                        do j = chebyshevDerivative%numCoeffs - 2, 1, -1
                                chebyshevDerivative%coeffs(j) = &
                                        chebyshevDerivative%coeffs(j + 2) &
                                       + 2.0_D * j * fit%coeffs(j + 1)
                        end do
                end if
                ll = changeX(chebyshevDerivative, chebyshevDerivative%ll)
                ul = changeX(chebyshevDerivative, chebyshevDerivative%ul)
                con = 2.0_D / (ul-ll)
                do j = 1, chebyshevDerivative%numCoeffs
                        chebyshevDerivative%coeffs(j) = &
                                chebyshevDerivative%coeffs(j) * con
                end do
                chebyshevDerivative%fitDone = .true.
        end function chebyshevDerivative
!C  (C) Copr. 1986-92 Numerical Recipes Software D041&0(9p#3.

        subroutine fitRationalFunction(rFit, func)

                use mConstants

                type(tRFit), intent(inout) :: rFit
                interface
                        function func(params, x)
                                use mPrecision
                                use mObject
                                real(D) :: func
                                type(tFitParams), intent(in) :: params
                                real(D), intent(in) :: x
                        end function func
                end interface
                real(D), parameter :: PIO2 = PI / 2.0_D
                integer :: numPts, numCoeffs, i
                real(D) :: ll, ul, hth
                real(D) :: realX
                real(D), target &
                   , dimension( NPFAC * (rfit%numOrder + rfit%denOrder + 1) ) &
                   :: xs, fs
                integer :: info

                numCoeffs = rFit%numOrder + rFit%denOrder + 1
                numPts = NPFAC * numCoeffs

                ll = changeX(rFit, rFit%ll)
                ul = changeX(rFit, rFit%ul)

                xGrid1: do i = 1, numPts / 2 - 1
!********************* This formula supposedly reduces roundoff errors at the
!********************* endpoints:
                        hth = PIO2 * (i - 1) / (numPts - 1)
                        xs(i) = ll + (ul - ll) * ( sin(hth) ) ** 2
                end do xGrid1
                xGrid2: do i = numPts / 2, numPts
                        hth = PIO2 * (numPts - i) / (numPts - 1) 
                        xs(i) = ul - (ul - ll) * sin(hth) ** 2
                end do xGrid2
                func_evals: do i = 1, numPts
                        if(rFit%xLog) then
                                realX = exp( xs(i) )
                        else
                                realX = xs(i)
                        end if
#ifdef VERBOSE
                        write(0, *) 'Calculating func(', realX, ')...'
#endif ! VERBOSE
                        fs(i) = func(rFit%params, realX)
#ifdef VERBOSE
                        write(0, *) fs(i)
#endif ! VERBOSE
                        if(rFit%yLog) fs(i) = log( fs(i) )
                end do func_evals

                rFit%coeffs = doFitRationalFunction(rFit, xs, fs)

#ifdef VERBOSE
                write(0, *) 'coeffs returned:'
                write(0, *) rFit%coeffs
#endif ! VERBOSE

                if(rFit%inverse) &
                        rFit%invCoeffs = doFitRationalFunction(rFit, fs, xs)

#ifdef VERBOSE
                write(0, *) 'inv returned:'
                write(0, *) rFit%invCoeffs
#endif ! VERBOSE

                rFit%fitDone = .true.

        end subroutine fitRationalFunction

        ! Can't have pure because of LA_GELS
        function doFitRationalFunction(rFit, x, y)

!                use SUNPERF ! Sun Performance Library (for subroutine DGELS)
                use LA_PRECISION, only: WP => DP
                use F95_LAPACK, only: LA_GELS

                type(tRFit), intent(in) :: rFit
                real(D), dimension(rfit%numOrder + rfit%denOrder + 1) :: &
                        doFitRationalFunction
                real(D), intent(in) :: x(:), y(:)

                ! Max # of least squares fits (5):
#ifdef DEV
                integer, parameter :: MAXIT = 5
#else ! DEV
                integer, parameter :: MAXIT = 50
#endif ! DEV

                real(D), dimension( size(x), 1 ), target :: &
                        coeffsData1, coeffsData2
                real(D), pointer, dimension(:, :) :: coeffs, scratchCoeffs
                real(D), dimension( size(x) ), target :: &
                        approxData1, approxData2
                real(D), pointer, dimension(:) :: approx, scratchApprox
                real(D), dimension( size(x) ) :: ee, wt
                real(D), dimension &
                         ( size(x), rFit%numOrder + rFit%denOrder + 1 ) :: u
                integer :: numPts, numCoeffs, i, it, j
                real(D) :: e, pow, sum
                real(D) :: relativeError, worstRelativeError, &
                           maxRelativeDeviation
                integer :: info
#ifdef VERBOSE
                real(D) :: xVal, yVal, yValApprox
#endif ! VERBOSE
!******************************************************************************
! for Sun only:
!************* Size of work array for least squares.
!************* After call, work(1) returns the optimum LWORK
!                integer, parameter :: LWORK = 512
!                real(D), target :: work(LWORK)
!Sun             nullify(ptrWork)
!Sun             ptrWork => work
!******************************************************************************

                nullify(coeffs, scratchCoeffs, approx, scratchApprox)

                numCoeffs = rFit%numOrder + rFit%denOrder + 1
                numPts = size(x)

                coeffs => coeffsData1
                scratchCoeffs => coeffsData2
                approx => approxData1
                scratchApprox => approxData2

                maxRelativeDeviation = huge(1.0_D)
                wt = 1.0_D
                ee = 1.0_D
                e = 0.0_D
                guess: do it = 1, MAXIT
                        do i = 1, numPts
                                pow = wt(i)
                                scratchCoeffs(i, 1) = &
                                    pow * (y(i) + sign( e, ee(i) ))
                                do j = 1, rFit%numOrder + 1
                                        u(i, j) = pow
                                        pow = pow * x(i)
                                end do
                                pow = -scratchCoeffs(i, 1)
                                do j = rFit%numOrder + 2, numCoeffs
                                        pow = pow * x(i)
                                        u(i, j) = pow
                                end do
                        end do
!****************************************************************************
! GELS is the standard least squares routine in the Sun library.
!                        call GELS('N', numPts, numCoeffs, 1, u, numPts, &
!                                  scratchCoeffs, numPts, ptrWork, LWORK, info)
!                        worksize: if( work(1) > LWORK ) then
!           write(0,*) 'INFO: Performance may be improved by increasing'
!           write(0,*) '      LWORK in FitRationalFunction to at least', work(1)
!           write(0,*) '      at least', work(1)
!           write(0,*)
!                                end IF worksize
!****************************************************************************
! LA_GELS is a least squares routine in the F95 interface to LAPACK
                        call LA_GELS(u, scratchCoeffs(1 : numPts, 1), 'N', info)
#ifdef DEBUG
                        checkInfo: if(info < 0) then
               write(0, *) 'ERROR: ', info, '-th argument to DGELS had an'
               write(0, *) '        illegal value'
               stop
                        end if checkInfo
#endif
!****************************************************************************
                        relativeError = 0.0_D
                        worstRelativeError = 0.0_D
                        sum = 0.0_D
                        find_max_error: do j = 1, numPts
                                scratchApprox(j) = &
                                        doEvalRationalFunction( x(j), &
                                            rFit%numOrder, rFit%denOrder, &
                                            scratchCoeffs(:, 1) )
                                ee(j) = scratchApprox(j) - y(j)
                                wt(j) = abs( ee(j) )
                                sum = sum + wt(j)
                                log_rel_err: if(rFit%xLog) then
                                        relativeError = &
                                            abs(exp( ee(j) ) - 1.0_D)
                                else log_rel_err
                                        if( y(j) /= 0.0_D ) &
                                            relativeError = abs( ee(j) / y(j) )
                                end if log_rel_err
                                if(relativeError > worstRelativeError) &
                                        worstRelativeError = relativeError
                        end do find_max_error
                        e = sum / numPts
                        if_better: if(worstRelativeError &
                            <= maxRelativeDeviation) then
                                maxRelativeDeviation = worstRelativeError
                                call switchPointers(coeffs, scratchCoeffs)
                                call switchPointers(approx, scratchApprox)
                        end if if_better
!                        write (6, 10) it, worstRelativeError
!10                      format(1x, 'FitRationalFunction iteration = ', I2, &
!                                ' max relative deviation = ', 1PE10.3)
                end do guess

                doFitRationalFunction = coeffs(1 : numCoeffs, 1)
#ifdef VERBOSE
                write(0, *) 'Coeffs:'
                write(0, *) coeffs(1 : numCoeffs, 1)

                write(0, *) 'Fit complete. Max deviation:', maxRelativeDeviation

                do i = 1, numPts
                        if(rFit%xLog) then
                                xVal = exp( x(i) )
                        else
                                xVal = x(i)
                        end if

                        yVal = changeY( rFit, y(i) )
                        yValApprox = changeY( rFit, approx(i) )

                        write(0, '(3e11.3)') xVal, yVal, yValApprox
                end do
#endif ! VERBOSE

        end function doFitRationalFunction

        real(D) pure function evalRationalFunction(rFit, x)
                type(tRFit), intent(in) :: rFit
                real(D), intent(in) :: x
                real(D) :: fitX
                fitX = changeX(rFit, x)
                evalRationalFunction = doEvalRationalFunction(fitX, &
                        rFit%numOrder, rFit%denOrder, rFit%coeffs)
                evalRationalFunction = changeY(rFit, evalRationalFunction)
        end function evalRationalFunction

        real(D) pure function fitEvaluateInverse(rFit,x)
                type(tRFit), intent(in) :: rFit
                real(D), intent(in) :: x
                real(D) :: fitX
                fitX = changeX(rFit, x)
                fitEvaluateInverse = doEvalRationalFunction(fitX, &
                        rFit%numOrder, rFit%denOrder, rFit%invCoeffs)
                fitEvaluateInverse = changeY(rFit, fitEvaluateInverse)
        end function fitEvaluateInverse

        !*****************************************************************
        ! Calculates the rational function 
        ! (p_0+p_1*x+...+p_m*x**m)/(1+q_1*x+...+q_k*x**k). The
        ! coefficients are stored as coeffs=(p_0,p_1,...,p_m,q_1,...,q_k).
        !*****************************************************************
        real(D) pure function doEvalRationalFunction(x, m, k, coeffs)
                real(D), intent(in) :: x
                integer, intent(in) :: m, k
                real(D), intent(in) :: coeffs(m + k + 1)
                integer :: j
                real(D) :: sumd, sumn
                sumn = coeffs(m + 1)
                do j = m, 1, -1
                        sumn = sumn * x + coeffs(j)
                end do
                sumd = 0.0_D
                do j = m + k + 1, m + 2, -1
                        sumd = ( sumd + coeffs(j) ) * x
                end do
                doEvalRationalFunction = sumn / (1.0_D + sumd)
        end function doEvalRationalFunction

        subroutine switch1DPointers(a,b)
                real(D), pointer :: a(:), b(:)
                real(D), pointer :: temp(:)
                temp => a
                a => b
                b => temp
        end subroutine switch1DPointers

        subroutine switch2DPointers(a,b)
                real(D), pointer :: a(:,:), b(:,:)
                real(D), pointer :: temp(:,:)
                temp => a
                a => b
                b => temp
        end subroutine switch2DPointers

        real(D) pure function rFitChangeX(rFit, x)
                type(tRFit), intent(in) :: rFit
                real(D), intent(in) :: x
                if(rFit%xLog) then
                        rFitChangeX = log(x)
                else
                        rFitChangeX = x
                end if
        end function rFitChangeX

        real(D) pure function rFitChangeY(rFit, y)
                type(tRFit), intent(in) :: rFit
                real(D), intent(in) :: y
                IF(rFit%yLog) then
                        rFitChangeY = exp(y)
                else
                        rFitChangeY = y
                end if
        end function rFitChangeY

        real(D) pure function cFitChangeX(cFit, x)
                type(tCFit), intent(in) :: cFit
                real(D), intent(in) :: x
                IF(cFit%xLog) then
                        cFitChangeX = log(x)
                else
                        cFitChangeX = x
                end if
        end function cFitChangeX

        subroutine checkAlloc(info)
                integer, intent(in) :: info
                if(info /= 0) then
                        write(0, *) 'WARNING: Allocation in'
                        write(0, *) '         FitRationalFunction failed!'
                end if
        end subroutine checkAlloc

        subroutine checkDealloc(info)
                integer, intent(in) :: info
                if(info /= 0) then
                        write(0, *) 'WARNING: Deallocation in'
                        write(0, *) '         FitRationalFunction failed!'
                end if
        end subroutine checkDealloc

end module mFit
